module time_int

  use globals
 
  implicit none
 
contains
!----------------------------------------------------------------

subroutine compute_dt

  implicit none
  
  real(8)  :: coeff
  
  ! RK4 stability coefficient, don't mess with unless blow-ups occur
  ! should be stable for all particles/cases...
  ! It would probably be worth updating this to implicit for very small
  ! particles.
  if (brun) then 
     coeff = 0.2d0 
  else
     coeff = 0.5d0 
  endif
  
  ! stability criteria based on smallest cell size
  dtg = min(dr**2.0d0 * coeff * (rho*cp)/(intCond*rad),dr**2.0d0 * coeff * (rho*cp)/(cond))
  dtl = tau/20.0d0
 
  nrep = floor(time*(rep*10.0d0**6.0d0))
  
  ! determine time step so that pulse is sufficiently resolved
!  if (time .lt. 4.0d0*tau .and. dtg .gt. dtl) then
  if (time .lt. tau .and. dtg .gt. dtl) then
	dt = dtl
  elseif (rep .gt. 0.0d0                                  .and. &
         & time+dtl .gt. (nrep/(rep*10.0d0**6.0d0))       .and. &
!         & time .lt. (nrep/(rep*10.0d0**6.0d0)+4.0d0*tau) .and. &
         & time .lt. (nrep/(rep*10.0d0**6.0d0)+tau)       .and. &
         & dtg .gt. dtl)                                  then
        dt = dtl
  else
	dt = dtg
  endif

end subroutine compute_dt 
!----------------------------------------------------------------

subroutine melt 

  implicit none
  
  if ((Ta .ge. 1337.0d0) .and. (Emelt .le. Ereq)) then
     Ta=1337.0d0
     Emelt=Emelt+dtp*Vp*Geph*(Te-Ta)-dtp*Qw  ! check this!
     if (mflag1 .eq. 0) then
        write(*,'(a,d11.5,a)') " >>> Particle melting at time ",time," <<< "
	mflag1 = 1
     endif
  endif
  
  if (Emelt .ge. Ereq .and. mflag2 .eq. 0) then
     write(*,'(a,d11.5,a)') " >>> Particle completely melted at time ",time," <<<"
     mflag2 = 1
  endif

end subroutine melt
!----------------------------------------------------------------

subroutine boil

  implicit none
  real(8)  :: delta_intCond, frac, decay
  
  decay = 0.01d0
  frac = 0.95d0
  delta_intCond = intCondo * (1.0d0 - 1.0d0/wvratio)
  
  if ((Tw(2) .ge. frac*Tw_spin)) then
     ! to avoid numerical instability from instant conductance change
     ! this effectivley "smears" the effects of phase transition.
     ! This should not be included in the governing eq as it is a physical discontinuity! 
     intCond = intCondo - delta_intCond*(1.0d0-exp(-decay*(Tw(2)-frac*Tw_spin)))
     if (bflag .eq. 0) then
        write(*,'(a,d11.5,a)') " >>> Water boiling at time ",time," <<< "
	bflag = 1
     endif
  else
     intCond = intCondo 
  endif

end subroutine boil
!----------------------------------------------------------------

! <warp>
subroutine update
  
  implicit none
  
  Ce = 70.0d0*Te
  
  Ca = 109.5790d0+0.1280d0*Ta-3.40d-4*(Ta**2.0d0) + &
     & 5.24d-7*(Ta**3.0d0)-3.93d-10*(Ta**4.0d0)   + &
     & 1.17d-13*(Ta**5.0d0)

  Ca = Ca*19300.0d0

  !!! HARD-CODING Ca HERE !!!
  Ca = 126.0D0 
	
  nrep = floor(timep*rep*10.0d0**6.0d0)
  ! offset of pulse energy for multiple pulses
  if (rep .gt. 0.0d0) then
     if (timep+dtp .gt. (nrep/(rep*10.0d0**6.0d0)) .and. timep.gt.0.0) then
        toffset = nrep/(rep*10.0d0**6.0d0)
     endif
  endif
	
  !!! MAKE THIS A READ-IN OPTION !!!

  ! Gaussian pulse, necessary for non-linear electron emission process
  ! where tau must be fwhm and Epulse is total integrated pulse energy
  ! this is power input (W)
  !Et = Epulse * 2.0d0/tau * (log(2.0d0)/pi)**0.5d0 * exp(-4.0d0*log(2.0d0)*((timep-2.0d0*tau-toffset)/tau)**2.0d0)

  ! Tophat pulse
  Et = 0.0D0
  if (time.lt.tau) Et = Epulse/tau
  
  ! print *, time, timep, dt, toffset, Et
  
end subroutine update
!----------------------------------------------------------------

subroutine dTe(Te_in,Ta_in,Te_out)

  use emit

  implicit none
  real(8) :: Te_in, Ta_in, Te_out
  real(8) :: eout, rdetot !, aeout_calc
  integer :: ii, dummy

  rdetot = detot
  if (eoutrun .and. rdetot .gt. 0.0d0) then
     eout = rdetot * aeout_calc(dummy)/dtp
  else
     eout = 0.0d0
  endif
  	
  ! modified qui-tien equation for electron heating with emission considered	
  Te_out = 1.0d0/(Vp*Ce) * (Vp*Geph*(Ta_in-Te_in)+Et-eout)
  
end subroutine dTe
!----------------------------------------------------------------

subroutine dTa(Te_in,Ta_in,Ta_out)

  implicit none
  real(8) :: Te_in, Ta_in, Ta_out	

  Qw = As*intCond*(Ta_in-Tw(1))  ! heat out from particle (W)
  ! qui-tien equation for phonon heating
  Ta_out = 1.0d0/(Vp*Ca) * (Vp*Geph*(Te_in-Ta_in)-Qw)
  
end subroutine dTa
!----------------------------------------------------------------
subroutine dTw(Tw_in,Tw_out)

  use space
  use emit
  
  implicit none
  integer :: ii, dummy
  real(8) :: alpha, esource, wetot, wetot_old !, ewater_heat
  real(8) :: aeout !, aeout_calc
  real(8), dimension(NI) :: Tw_out, Tw_in

  ! spherical diffusion, no convection
  alpha = cond/(rho*cp) ! m^2/s
  Qw_p = As*intCond*(Ta_old-Tw_in(1))  ! heat out from particle (W)
  
  if (eoutrun) then
     aeout = aeout_calc(dummy) - bh*eVJ
  else
     aeout = 0.0d0
  endif  
  
  ii = 1 ! particle bc
  if (eoutrun) then
     wetot = etot
     esource = ewater_heat(ii, wetot, wetot_old, aeout)
     wetot_old = wetot
  else
     esource = 0.0d0
  endif
  
  call qdiff ! find Q(r)
  Tw_out(ii) = -1.0d0/(4.0d0*pi*r(ii)**2.0d0*rho*cp) * fdiff1(Qr,r,ii)
  
  do ii=2,NI-1 ! rest of domain
     if (eoutrun) then
	wetot = etot
        esource = ewater_heat(ii, wetot, wetot_old, aeout)
	wetot_old = wetot
     else
        esource = 0.0d0
     endif

     Tw_out(ii) = alpha/r(ii) * (2.0d0*diff1(Tw_in,r,ii) + &
	        & r(ii)*diff2(Tw_in,r,ii)) + esource/dt
  enddo
    
end subroutine dTw
!----------------------------------------------------------------
subroutine check_dTw(time_c,tstep_c)

  use write_results
  implicit none
  real(8) :: time_c
  integer :: tstep_c
  
  ! terminate simulation if heating of three systems is complete
!  if (time_c .gt. 4.0d0*tau .and. Tw(2) .le. Tw_last) then
  if (time_c .gt. tau .and. Tw(2) .lt. Tw_check) then
     write(*,*) " >>> Heating complete at time: ", time_c, " <<< "
     call write_T(tstep_c,1)
     stop
  endif
  
end subroutine check_dTw
!----------------------------------------------------------------

subroutine particle_step

  use emit

  implicit none
  integer :: subt
  real(8) :: rsubt, rtpstep
  real(8) :: Te_temp, Ta_temp
  real(8) :: Te1, Te2, Te3, Te4
  real(8) :: Ta1, Ta2, Ta3, Ta4
  
  ! For stable substep for particle heating, somewhat arbitrary 
  ! for now (conservatively).  Should be based on either pulse 
  ! duration or e-p coupling constant.
  subt = 10
  rsubt = subt
  dtp=dt/rsubt 
  timep = time
  Ta_old = Ta
  Te_last = Te
  Ta_last = Ta
  
  do tpstep=1,subt
	
	rtpstep = tpstep
	timep = timep + dtp
	call update                 ! update Ce, Ca, and Et
	E_total = E_total + Et*dtp  ! total laser energy, just a sanity check

	heatLoss = heatLoss+dtp*Qw
	heatIn = heatIn+dtp*Et
       
	! begin RK4
	call dTe(Te,Ta,Te1)
	call dTa(Te,Ta,Ta1)
	   
	Te_temp = Te + Te1*0.5d0*dtp
	Ta_temp = Ta + Ta1*0.5d0*dtp
	call dTe(Te_temp,Ta_temp,Te2)
	call dTa(Te_temp,Ta_temp,Ta2)
	   
	Te_temp = Te + Te2*0.5d0*dtp
	Ta_temp = Ta + Ta2*0.5d0*dtp
	call dTe(Te_temp,Ta_temp,Te3)
	call dTa(Te_temp,Ta_temp,Te3)
	   
	Te_temp = Te + Te3*0.5d0*dtp
	Ta_temp = Ta + Ta3*0.5d0*dtp
	call dTe(Te_temp,Ta_temp,Te4)
	call dTa(Te_temp,Ta_temp,Te4)
	   
	Te = Te + dtp*1.0d0/6.0d0*(Te1+2.0d0*Te2+2.0d0*Te3+Te4)
	Ta = Ta + dtp*1.0d0/6.0d0*(Ta1+2.0d0*Ta2+2.0d0*Ta3+Ta4)

    ! if you wish to consider phase transition of particle
    if (mrun) call melt 
	
    ! calculate Efield and electron emission
    if (erun) call emission

    Ta_last = max(Ta_last,Ta)
    Te_last = max(Te_last,Te)
	   
  enddo
  
end subroutine particle_step
!----------------------------------------------------------------

subroutine water_step
  
  use space
  implicit none
  integer  :: ii
  real(8), dimension(NI) :: Tw_temp
  real(8), dimension(NI) :: Tw1, Tw2, Tw3, Tw4

  Tw_check = Tw(2)
  Tw_last = Tw(1)
  Tw(NI) = Tamb  
  Tw_temp(NI) = Tamb 
  !Qw_p = As*intCond*(Ta_old-Tw(1))  ! heat out from particle (W)

  ! begin RK4
  call dTw(Tw,Tw1)
  do ii=1, NI-1
	Tw_temp(ii) = Tw(ii) + Tw1(ii)*0.5d0*dt
  enddo
  
  call dTw(Tw_temp,Tw2)
  do ii=1, NI-1
	Tw_temp(ii) = Tw(ii) + Tw2(ii)*0.5d0*dt
  enddo
  
  call dTw(Tw_temp,Tw3)
  do ii=1, NI-1
	Tw_temp(ii) = Tw(ii) + Tw3(ii)*0.5d0*dt
  enddo
  
  call dTw(Tw_temp,Tw4)
  do ii=1, NI-1 
	Tw(ii) = Tw(ii) + dt*1.0d0/6.0d0*(Tw1(ii)+2.0d0*Tw2(ii)+2.0d0*Tw3(ii)+Tw4(ii))
  enddo

  Tw_last = max(Tw_last,Tw(1))
  Tw_check = max(Tw_check,Tw(2))
  
  ! bubble stuff would be added here, see lattice-boltzmann or level-set methods...
  if (brun) call boil

end subroutine water_step
!----------------------------------------------------------------

subroutine time_step

  use write_results
  implicit none
  
  real(8) :: rtstep, rtpstep
  
  ! intialize a bunch of stuff
  intCond = intCondo
  heatLoss = 0.0d0
  heatIn = 0.0d0
  Emelt = 0.0d0
  time = 0.0d0
  tstep = 0
  Ta = Tini
  Te = Tini
  Tw = Tini
  
  mflag1 = 0
  mflag2 = 0
  etot = 0
  eflag1 = 0
  eflag2 = 0
  eflag3 = 0
  bflag = 0
  nrep = 0
  toffset = 0.0d0
  Tw_flag = 0
  Ta_flag = 0
  Te_flag = 0

  !do tstep=0,maxstp      ! begin loop over timesteps
  do while (time .lt. tmax)
    tstep = tstep + 1
  
    call compute_dt      ! find stable & well-resolved step
    rtstep = tstep       ! these are just int to real conversions
	
    call particle_step   ! calculate particle temp
    call water_step      ! calculate water temp
    time = time + dt     ! update time
	
    call write_T(tstep,0)! write results
    if (efficient)  call check_dTw(time,tstep) ! avoid wasted computation if heating complete
			  
  enddo
  
end subroutine time_step 
!----------------------------------------------------------------
 
end module time_int
